Towards absorbing outer boundaries in General Relativity 



Luisa T. Buchman^^'^ and Olivier C. A. Sarbach^^'^'* 
^ Relativistic Astrophysics 169-327, Jet Propulsion Laboratory, 
4800 Oak Grove Drive, Pasadena, California 91109-8099, USA 
^Theoretical Astrophysics 130-33, California Institute of Technology, 
1200 East California Boulevard, Pasadena, California 91125-0001, USA 
'^Department of Mathematics, University of California at San Diego, 
9500 Gilman Drive, La Jolla, California 92093-0112, USA and 
'^Instituto de Fisica y Matemdticas, Universidad Michoacana de San Nicolas de Hidalgo, 
Edificio C-3, Cd. Universitaria. C. P. 58040 Morelia, Michoacdn, Mexico 
(Dated: February 7, 2008) 

We construct exact solutions to the Bianchi equations on a flat spacetime background. When the 
constraints are satisfied, these solutions represent in- and outgoing linearized gravitational radiation. 
We then consider the Bianchi equations on a subset of flat spacetime of the form [0, T] x Bji, where 
Br is a ball of radius R, and analyze different kinds of boundary conditions on dBa. Our main 
results are: i) We give an explicit analytic example showing that boundary conditions obtained 
from freezing the incoming characteristic flelds to their initial values are not compatible with the 
constraints, ii) With the help of the exact solutions constructed, we determine the amount of 
artiflcial reflection of gravitational radiation from constraint-preserving boundary conditions which 
freeze the Weyl scalar '^q to its initial value. For monochromatic radiation with wave number k and 
arbitrary angular momentum number I > 2, the amount of reflection decays as {kR)~'^ for large 
kR. iii) For each L > 2, we construct new local constraint-preserving boundary conditions which 
perfectly absorb linearized radiation with £ < L. (iv) We generalize our analysis to a weakly curved 
background of mass M, and compute first order corrections in M/R to the reflection coefficients 
for quadrupolar odd-parity radiation. For our new boundary condition with L = 2, the reflection 
coefficient is smaller than the one for the freezing ^Pq boundary condition by a factor of M/R for 
kR > 1.04. Implications of these results for numerical simulations of binary black holes on finite 
domains are discussed. 

PACS numbers: 04.20.-q, 04.25.-g, 04.25. Dm 



I. INTRODUCTION 



A common approach for numerically solving the Einstein field equations on a spatially unbounded domain is to 
truncate the domain via an artificial boundary, thus forming a finite computational domain 17 with outer boundary 
In order to obtain a unique Cauchy evolution, it is necessary to impose boundary conditions at d^l. These 
boundary conditions should form a well posed initial boundary value problem (IBVP) and, ideally, be completely 
transparent to the physical problem on the unbounded domain. Short of achieving the ideal, one can try to develop 
so-called absorbing boundary conditions which form a well posed IBVP and insure that only a very small amount 
of spurious gravitational radiation is reflected from dU, into the computational domain. Once the IBVP on f7 is 
formulated, it is solved via a numerical approximation scheme which, together with the truncation of the domain, 
introduces two artificial parameters: a discretization parameter h, describing the coarseness of the discretization, and 
a cut-off parameter R, which gives the size of the spatial domain $1. For a stable discretization, it is expected that 
the continuum solution of the unbounded problem is recovered in the limit where h ^ Q and R ^ oo. In practice, 
due to finite computer resources, it is not possible to take this limit. Instead, one needs to quantify how small h and 
how large R need to be so that the error is below a certain tolerance value. 

In this article, we address the "i?-dependent" part of this task. We analyze boundary conditions which have been 
recently presented in the literature, and provide estimates for the amount of spurious radiation coming from dfl. 
Additionally, we propose new boundary conditions for Einstein's vacuum field equations which introduce significantly 
less reflections than existing conditions. 

There has been a substantial amount of work on the construction of absorbing (also called non-refiecting in the 
literature) boundary conditions for wave problems in acoustics, electromagnetism, meteorology, and solid geophysics 
(see 0] for a review). One approach is based on a sequence of local boundary conditions with increasing order 

of accuracy. Although higher order local boundary conditions usually involve solving a high order differential equation 
at the boundary, the problem can be dealt with by introducing auxiliary variables at the boundary surface [5lja]. A 
different approach is based on fast converging series expansions of exact nonlocal boundary conditions (see Q and 
references therein). Of particular interest for this article is the work by Lau [MSEIj which generalizes the work in 
Ref. to the construction of exact non-reflecting boundary conditions for the Regge- Wheeler and Zerilli equations, 
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describing linear gravitational fluctuations about a Schwarzschild black hole. This approach is robust, very accurate, 
and stable. However, it is based on a detailed knowledge of the solutions which might not always be available in more 
general situations. 

For the fully nonlinear Einstein equations, the construction of absorbing outer boundary conditions is particularly 
difficult. First of all, Einstein's field equations determine the evolution of the metric tensor, so one does not know 
the geometrical structure of the spacctime before actually solving the IBVP. Hence, it is not clear a priori how the 
geometry of the outer boundary evolves. This poses a problem if one wants to fix, for example, the area of the boundary 
(917 to its initial value. Second, in the Cauchy formulation of Einstein's field equations, there exist constraint-violating 
modes which propagate with nontrivial characteristic speeds. This is in contrast to the standard Cauchy formulation 
of Maxwell's equations, where the evolution equations imply that the constraint variables (namely, the divergence of 
the electric and magnetic fields) are constant in time. Since the constraint variables in General Relativity propagate 
non-trivially, constraint-preserving boundary conditions (CPBC) must be specified so that constraint violations are 
not introduced into the computational domain. Finally, in General Relativity, it is difficult to define precisely what 
is meant by outgoing and ingoing radiation. This is due to the nonlinear nature of the theory and its diffcomorphism 
invariance. (See Ref. [Tlj for a discussion of this problem for nonlinear gravitational plane waves.) These issues all 
contribute to the challenge of determining the amount of spurious reflections from the outer boundary. 

A signiflcant advance towards developing absorbing boundary conditions for General Relativity was the first (and, 
to date, the only) well posed IBVP for Einstein's vacuum field equations presented in Ref. M^. This work, which 
is based on a tetrad formulation, recasts the evolution equations into a first order symmetric hyperbolic form with 
maximally dissipative boundary conditions, for which (local in time) well posedness is guaranteed [iJI • The boundary 
conditions constructed in [T^ control part of the geometry of the boundary surface by specifying its constant mean 
curvature, control the radiation by prescribing suitable combinations of the complex Newman-Penrose scalars a-nd 
\l/4, where the null tetrad is constructed from the evolution vector field and the normal to the boundary, and are 
constraint-preserving. Recently, there has been considerable effort to generalize the work in Ref. Il2(| by specifying 
CPBC for the more commonly used metric formulations of gravity (see Refs. mill EE m mili & eh and 
references therein). In particular, the methods in Refs. [igI \vh. ITM Il9l l2]| in addition to preserving the constraints, 
regulate the dynamical degrees of freedom by freezing the Newman-Penrose scalar defined with respect to a 
suitably chosen null tetrad, to its initial value 77] . 

Work focused on eliminating reflections from the outer boundary during fully relativistic vacuum simulations has 
been performed by several authors. In Ref. |22| . boundary conditions based on the work in Ref. which are perfectly 
absorbing for quadrupolar solutions of the flat wave equation, are numerically implemented via spectral methods, and 
used in a constrained evolution scheme of Einstein's field equations |2^. In Refs. 011^, solutions of the full nonlinear 
Einstein equations on a finite computational domain are matched to exact analytic, purely outgoing solutions of the 
weak field equations at the outer boundary of the domain. Refs. [26l WK . I2M |29j generalize this idea by matching 
the nonlinear equations to an "outer module", a code in which the equations are linearized about a Schwarzschild 
background, in order to carr y th e waveforms far into the wave zone. However, at the interface where the matching 
occurs, the methods in Refs. |24. 125ll2a . I27ll28| do not take into account either the constraints of the nonlinear Cauchy 
code or the characteristic structure of the nonlinear evolution equations, so it is not clear if the resulting problem is 
well posed (the work in 29], on the other hand, does take into account the constraints and the characteristic fields of 
the Cauchy code, and gives an implementation for the spherically symmetric Einstein equations coupled to a massless 
scalar field). Two other approaches presented in the literature for constructing absorbing boundary conditions are: 
matching the nonlinear Cauchy code to a nonlinear characteristic code (see Ref. for a review and ^3 for recent 
work) and matching an incoming characteristic formulation to an outgoing one at a time-like cylinder l33l . Finally, 
methods that avoid introducing an artificial outer boundary altogether compactify spatial infinity [s^. l34l]. or make 
use of hyperboloidal slices and compactify null infinity (see, for instance, 35, 36, 37]). 

In this article, we take a step closer to the construction of absorbing boundary conditions in General Relativity. In 
order to do so, the IBVP of Einstein's field equations is analyzed on a compact domain J7 C with smooth outer 
boundary dVL and two simplifying assumptions. The first assumption is that at all times, the boundary surface dVL is 
far from the strong field region, so that the gravitational field near the outer boundary is weak. As a consequence, the 
field equations can be linearized to a first approximation about fiat spacetime in the vicinity of the outer boundary. 
The linearized field equations can be conveniently described by the Bianchi equations, which yield a Lorcntz-invariant 
system for the linearized Weyl tensor having a structure which is very similar to that of Maxwell's equations. Moreover, 
the linearized Weyl tensor is invariant with respect to infinitesimal coordinate transformations, since it vanishes on 
the background |38| . so there are no gauge modes. The second assumption is that the boundary dVl. is approximately 
a metric sphere of area AttR'^. This assumption is quite natural. In fact, modern numerical relativity codes based on 
multi-block finite differencing js^ l40l l4lj ] or pseudo-spectral methods p3 . are designed to handle spherical outer 
boundaries. 

Under these assumptions, it is sufficient to analyze the Bianchi equations on a domain = Bji consisting of a ball 
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of radius R. We can then conveniently expand the hnearized Weyl tensor in terms of spherical tensor harmonics, 
because of the spherical symmetry of Bj^. The resulting equations are decoupled: they are a family of partial 
differential equations in one spatial dimension parameterized by the angular momentum numbers £ and to. For each 
fixed £ and m, the purely dynamical degrees of freedom can be described by a master equation for the Ncwman- 
Penrose scalar This equation admits exact solutions which propagate along either in- or outgoing null radial 
geodesies. The in- and outgoing solutions are related to each other by a time reversal symmetry t i—f ~t, making it 
possible to define sensibly in- and outgoing gravitational radiation. Hence, in our setting, it is clear how to quantify 
the amount of spurious radiation reflected at dBn. Using these exact solutions, we analyze the quality of boundary 
conditions which have been proposed in the literature; namely, those which freeze all the incoming characteristic fields 
to their initial values, and CPBC which freeze the Weyl scalar '^q to its initial value. Furthermore, we offer a set 
of improved CPBC, which are perfectly absorbing for linearized radiation on a Minkowski background up to some 
arbitrary multipolc number £. Finally, we extend our analysis to a weakly curved background. 

Our main results are the following. First, we show that the naive boundary condition which freezes all the incoming 
characteristic fields to their initial values is not compatible with the constraints. To show this, we construct explicit 
solutions to the IB VP which have the property that they satisfy the constraints exactly on the initial time slice t — 0, 
but violate them at later times t > 0. Second, we impose CPBC and freeze the Weyl scalar ^Pq at the boundary to its 
initial value. The exact outgoing solutions do not satisfy this boundary condition exactly. Specifically, the quantity 
^0 constructed from these solutions falls off as along the null geodesies t — r + const., where r denotes the areal 
radius coordinate. This means that a solution to the IBVP corresponding to the boundary condition 9t^'o — consists 
of a superposition of an in- and an outgoing wave, where the magnitude of the ratio of the ingoing to the outgoing 
wave amplitudes (which we define as the reflection coefficient) measures the amount of spurious reflection. We find 
that for monochromatic radiation with wave number k and arbitrary angular momentum number £ > 2, the reflection 
coefficient decays as {kR)~'^ for large kR. In particular, the reflection coefficient lies below 0.1% for quadrupolar 
radiation with kR > 6.4. Third, for each L > 1, we construct local CPBC Bl which, for L > 2, improve the CPBC 
involving 9t^'o = 0, being perfectly absorbing for linearized gravitational radiation on Minkowski space with angular 
momentum number £ < L. (Bi is just the freezing ^'o boundary condition and there is no improvement.) Since in 
many practical situations one expects the few lower multipolcs to dominate, an implementation of Bl for L = 2, 3, or 
4 should result in only a small amount of spurious reflection. For L = 2, our improved boundary condition B2 reads 

dt{dt + dr){r''^o)l^j, = 0. (1) 

Finally, we take into account flrst order corrections from the curvature of the background. Since we assume that 
the outer boundary lies in the weak field regime, we describe spacetime near the outer boundary by a perturbed 
Schwarzschild metric of mass M, thereby generalizing our previous analysis by taking into account curvature near the 
outer boundary. To estimate the effects due to curvature, we compute the flrst order corrections in 2M/R to the exact 
in- and outgoing solutions with £ = 2 and odd parity, and then recalculate our reflection coefficient for the CPBC 
involving dt'^o = 0. We flnd that for 2M/R <C 1, the corrected £ = 2 odd-parity reflection coefficient depends only 
weakly on 2M/R. In fact, our results indicate that the reflection coefficient even decreases when 2M/R increases (but 
stays small). For quadrupolar solutions satisfying the improved boundary condition (Q, which is perfectly absorbing 
for M = 0, we find that the reflection coefficient decays as {2M / R){kR)^^ for large kR and small 2M/R. More 
precisely, the reflection coefficient is less than the one for the freezing ^Pq boundary condition by a factor of M /R for 
kR > 1.04. 

This work is organized as follows. In Sect. ^ we write down the Bianchi equations on an arbitrary spacetime. 
These equations can be obtained from the Bianchi identities after imposing the Einstein field equations. Next, we 
assume the existence of a spacelike foliation and a preferred radial direction in each timeslice and perform a 2 -I- 1 -I- 1 
split of the Bianchi equations, which separate into evolution and constraint equations. The constraint propagation 
system describing the evolution of constraint errors is also discussed. 

In Sect, mil we specialize to a flat spacetime background of the form [0, T] x S^, where Br denotes a ball of 
radius R. By performing a decomposition into spherical tensor harmonics with angular momentum number £, we 
show that for each £ > 2, the Bianchi equations can be reduced to two master equations. The flrst master equation 
describes the propagation of constraint violations, and is homogeneous. The second is an equation for 5*2, describing 
the propagation of linearized gravitational radiation, and has a source term which depends on the solution of the first 
master equation. One of the advantages of working with a master equation for ^'2 instead of a master equation for ^'o 
or ^"4, as is usually done when studying perturbations of black holes with a Petrov type D metric 43, 44J, is that for 
linearization about Minkowski spacetime, the former is invariant with respect to time reversal. Consequently, there 
is a nice symmetry between in- and outgoing solutions: one can be obtained from the other by changing the sign 
of t. This symmetry makes it possible to define the refiection coefficients in a natural way. In contrast, under time 
reversal, is mapped to conjugate ^'4 and vice versa, so that the in- and outgoing parts of ^'o look quite different. It 
is shown in this section that the master equations governing the constraint violations and the gravitational radiation 
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both admit exact analytical solutions, which can be obtained by applying suitable differential operators to the solution 
of the one-dimensional flat wave equation. These solutions can be split in a unique way into in- and outgoing solutions 
describing, respectively, in- and outgoing constraint violations or in- and outgoing gravitational radiation. 

In Sect. IIVI we use the exact in- and outgoing solutions found in the previous section to construct exact solutions 
to the IB VP on Bn corresponding to different boundary conditions on dBj^. In Sect. IIV Al we start by analyzing the 
characteristic structure of the evolution equations, and specify boundary conditions which freeze the incoming fields 
to their initial values. The incoming fields are related to the Weyl scalars and ^'i, so these boundary conditions 
freeze 4'o and 5*1 at the outer boundary to their initial values. By constructing an explicit solution with constraint 
satisfying data at i = 0, we show that these "freezing" boundary conditions are not compatible with the constraints 
in the sense that the solution violates the constraints for t > 0. Next, in Sect. IIV Bl we replace the boundary 
condition which freezes ^'i to its initial value with CPBC which guarantee that solutions of the IBVP satisfy the 
constraints everywhere on Bn and at all times t > 0, provided they hold initially. This can be achieved in two 
ways. The first is the one proposed in Ref. |l2j| . which adds suitable combinations of the constraint equations to the 
evolution equations so that at the boundary, the constraints propagate tangentially to the boundary. The second is 
to analyze the characteristic structure of the constraint propagation system, and set the incoming constraint fields 
to zero at the boundary. Assuming in what follows that the constraints are satisfied exactly, we consider only the 
homogeneous master equation for '^2- We impose the freezing boundary condition dt'^o — at the boundary r — R 
on a superposition of in- and outgoing monochromatic waves for arbitrary £, and calculate the resulting reflection 
coefficients. These coefficients, which depend only on the dimensionless quantity kR (where k is the wave number) 
are of order unity if kR < i, and decay as {kR)~^ for large kR. In Sect. IIV CI we construct the hierarchy Bi, 82,- ■■ 
of improved boundary conditions, having the property that Bl is perfectly absorbing for all linearized gravitational 
waves with angular momentum number € up to and including L. The construction of these boundary conditions is 
strongly related to the hierarchy proposed in 

Finally, in Sect. ^ we generalize our analysis to odd-parity perturbations of a Schwarzschild background of mass 
M . Assuming that M/R <C 1, we compute first order corrections in M/R to the reflection coefficient corresponding 
to the freezing boundary condition, for 1 = 2. In addition, we compute the reflection coefficient for the boundary 
condition B2 (which is perfectly absorbing for M = 0), and show that it is smaller than the one for the freezing 
condition by a factor of M/R for kR > 1.04. 

Implications for the modeling of isolated systems such as a binary black holes are discussed in the conclusions. In 
an appendix, we show that the IBVP corresponding to the master equation for '^2 and our new boundary conditions 
S2, B3,... is stable in the sense that the solutions depend uniquely and continuously on the initial data. 



II. THE BIANCHI IDENTITIES 



We consider the Bianchi equations 

VaC^bcd — Jbcd (2) 

on a given background geometry {M,gab), with Cabcd a tensor field possessing the same algebraic symmetries as the 
Weyl tensor: 

C[abc]d = 0, C[ab][cd] — Cabcd = Ccdab , Q^'^Cabcd = 0, (3) 

and Jbcd a given source tensor which is traceless and satisfies J[bcd] — Jb[cd] = J>.^w|78j. Eq. ^ has its origin in the 
Bianchi identities, 

^aW\cd = V[c (Gd]b - lgd]b9''^Gef^ , (4) 

where Wabcd and Gab denote, respectively, the Weyl tensor and the Einstein tensor belonging to the metric gab- If 
Einstein's equations are imposed, then the right-hand side of the identity (0J can be re-expressed in terms of the 
stress-energy tensor, and Q becomes an equation for the Weyl tensor which is of the form of Eq. The identity 

VV7 /^ab /-tab tdc 

aVfeO cd — (-/ e[c^ d]ab , 

where Rabcd denotes the Riemann tensor belonging to the background metric gab, yields the integrability condition 

^^Jbcd = C"^'^ c\cR^ d]ab ■ (5) 
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The right-hand side of this equation vanishes if gab is flat or conformally flat. 

In Sections lllll and HVI we will assume that the background geometry {M,gab) is flat, in which case Eq. Q 
describes the propagation of linearized gravitational radiation, with Cabcd the linearized Weyl tensor. Since the Wey\ 
tensor vanishes for flat spacetime, Cabcd is invariant with respect to infinitesimal coordinate transformations [38j. 
As a consequence, the Bianchi equations are well-suited for studying linearized gravitational waves since they are 
manifestly gauge-invariant. For a flat background geometry, the integrability condition Q reduces to the requirement 
that Jbcd be divergence-free. 

Finally, the Bianchi equations Q can be coupled either to equations for metric components and Christoffel symbols, 
or to equations for tetrad fields and connection coefhcients, giving the full nonlinear vacuum Einstein equations 

[n mil mm 

A. 3-1-1 split 

We assume there exists a globally defined time function i : M ^ R such that M is foliated by spacelike hypersurfaces 
'Et = {p € M : t{p) — t}. Let — — aVai be the future-directed unit normal to these slices, where the time 
orientation is chosen so that the lapse function, a, is strictly positive. The three-metric hab and extrinsic curvature 
kab are defined asjT^ 

hab = gab + naTib , kab = Va?T.6 + nattb , 

where at, = n°'V anb is the acceleration along the integral curves of n^. For a one-form Va tangential to E,- in the sense 
that ^'af^° — 0, the spatial covariant derivative DaVb is defined as DaVb = ha'^hb'^W cVd- The spatial covariant derivative 
of a general tangential tensor field is defined similarly. The electric and magnetic parts of Cabcd are, respectively. 



Eab — Cacbdn'^n'^^ Hab — Ti^'^^caefS 



2 - ' 

where Zbcd = n'^Sabcd denotes the natural volume element on (S^, hab)- From the symmetries © of Cabcd, it follows 
that Eab and Hab are symmetric, traceless, and orthogonal to n". Furthermore, the ten fields {Eab, Hab} uniquely 
determine Cabcd- 

Cabcd — — 4n[ai?b][crici] — Sab^EefE-^ cd — '^n^aHb\e£'^ cd + 2£afc'^-ffe[c"-d] • 

The decomposition of Eq. Q into components normal and tangential to yields the evolution equations 

£nEab - ~e,,d(a{D^ + 'ia')H\) +bk(^a''Ei,)d'2kEab - habk"'E,d + Rab , (6) 
£nHab = +e,dia{D' + 2a')E\j+5k(^a''Hi,)d-2kHab-habk'''H,d + Sab, (7) 



and the constraint equations 



In these equations, k = h ka, 



D'Eab - k'"'s\aH,b = Pa , (8) 
D^'Hab + k'^'e^aEcb = Qa ■ (9) 



Pc = n'ri'^Jbcd , Qa = -l: n^'Sa^'^Jbcd , 



2 

Rcf = —hf^^rfhf-^''' Jbcd , Sef = — - hf^ J' £ j-^"^"^ Jbcd , 



and £n denotes the Lie derivative with respect to the unit nor mal field 7i°. Notice that Eqs. (P7|) and obey the 
"Dirac duality" symmetry 

[Eab, Hab) ^ {Hab, -Eab), {Pa, Qa) {Qa, -Pa), {Rab, Sab) >-> {Sab, - Rab) - (10) 

B. 2-1-1 split 

In addition to the foliation by spacelike hypersurfaces, the existence of a unit spatial vector field s° which is 
everywhere tangential to the hypersurfaces St is assumed. The existence of such a vector field allows us to introduce 
a Newman-Penrose null tetrad 

r^^K + s"^), k" ^ ^ {n" - s") , to" = ^K-t-iU."^), = ^{v" -iw"), 

y/2 V2 V2 V2 
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where w° and w°' are two mutually orthogonal unit vector fields which are normal to and s". The corresponding 
Newman- Penrose Weyl scalars |43 are defined as jSOj 

Next, we decompose Eab and iJaj, into components parallel and normal to Sa- More precisely, we write 

1 



Eab 



SaSb - -jlab ) E + 2S(aEb) + Eab , 



1 



Hab = [ SaSb - -Jab ) H + 2s(ai/fc) + Hab 



where jab = hab - SaSb, E = EabS°-s^, Ea = ja^'Ebcs", and Eab = {la^Jb" 
H, Ha and Hab)- In terms of these quantities, the Weyl scalars are 



*0 
^'4 



Eab + ea^Hcb mrm" 

[Ea+Sa'Hb] ™^ 

l[E-^H], 
-^[Ea-ea'Hb] rh", 



Ea 



Sa'H, 



cb 



Ecd (with similar expressions for 

(11) 
(12) 

(13) 
(14) 

(15) 



where Eab = Sabcs"- 

In order to decompose the evolution equations (|6I7|I . we make additional assumptions on the vector fields n° and 
s". First, we assume that s" is geodetic and everywhere orthogonal to closed 2-surfaces Sr in Sj. This implies that 

DaSb — f^ab i 

where Kab is a symmetric tensor field which is orthogonal to s", representing the extrinsic curvature of the 2-surfaces 
Sr as embedded in Et. Next, we assume that the Lie-derivative £ns'^ of the vector field s° with respect to can be 
written as a linear combination of n° and s". This implies that any covariant tensor field ta^a-z.-.a^ which is orthogonal 
to n" and s" has the property that £ntaia2...a^ is again orthogonal to and s°. Finally, we assume that the Lie- 
derivative £nSa of the one-form Sa is proportional to Sa- These properties, together with the relations Jias" — 0, 
Sas" = 1, £nna = Da{loga), and £nhab = 2A;afc, imply that 



la hcS'' 



0, 



^ns" = (£s loga)n° - A:s°, £„Sa = ksa 



where k = kabS°'s^ . Although these assumptions are strong, and may not all hold for a generic spacetime, they 
are satisfied for the background spacetimes and foliations used in this article. In particular, they arc met for any 
spherically symmetric spacetime of the form 



^dt^ + 7^ {dr + 13 dtf + (di?^ + sin^ d dip^) , 



where a, (3 and 7 are smooth functions of t and r, and where Uadx"^ = —adt and Sadx"^ — j{dr + (3dt). Decomposing 
the extrinsic curvature kab as 



kab 



1 



1 



SaSb - Tilab ]k + kab + Tilabk 
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where kab = {'Ja'^Jb'^ — 57ah7'^'') kcd, we find from evolution equation © that 



£nE = 2 ^"'^^ i^^^b) + k^'Sa^Hcb +-{k-k)E~ k'^^'Eab + R , 



2a2 



2a2 



4a2 



2kJ'eb'Hc + \{k~ 3fc) -Ba + + Ra , 

- (5A: - fc) i^ab - ^fcab^^ + 5k(^a''Eb)c ~ iflabk^'^Ecd + Rab , 



£{a''fib)cH 



tf 



(16) 



(17) 



(18) 



where [■■■Y-^ denotes the trace-free part with respect to jab, >^ab} denote the trace and trace-free part of Kab, 
respectively, V denotes the covariant derivative compatible with jab, and {R, Ra, Rab} denote the parallel/parallel, 
parallel/transverse, and transverse trace-free parts of Rab, respectively. The constraint equation © yields 



P ^ £sE + V'Ea - k'^'Eab + + k^-ha'Hbc , 

Pa = £sEa + V'Eab - \VaE + KEa ~ ^ {Sk - k) Ea'^Hb " ka'^Sb'H, 



(19) 
(20) 



where we have defined P — PaS°' and Pa — ja^Pb- Evolution and constraint equations for H, Ha and Hab are easily 
obtained by applying the Dirac duality transformations (|10|1 . 



C. Propagation of the constraint fields 

The decomposition of the integrability condition ((SJ into parts normal and tangential to gives (assuming that 
the background metric is flat or conformally flat) 

£nPa = -\ea''\D, + 3ac)Qd + I {ka'Pb - kPa) + {D' + a'')Rab + ea'"'k\Sbd , (21) 

£aQa = \ea'''{D,+?,a,)Pd + ^{ka''Qb~kQa)+{D'' +a'')Sab~ea'''k\Rbd. (22) 

If the evolution equations H6I7|I hold with Rab — Sab — 0, then the fields Pa and Qa obey homogeneous Maxwell-like 
equations with transmission speed half the speed of light. In particular, Pa = Qa = on an initial spatial slice gives 
-Pa = Qa = on the future domain of dependence of the initial slice. We may therefore regard Pa = Qa = as 
constraints which are propagated by the evolution equations (|6l7ll with Rab = Sab = 0. 

Under the same assumptions on the vector fields and as in the previous subsection, the 2-1-1 split of the 
constraint propagation equations yields 

£nP = -^e''''Va{a^Qb) + l{k-3k)P 

+ -£s{aR) + -V\aRb)~k''''Rab + ^R + k''''ea'Sbc, (23) 

£nPa = ^ea'[£s{a^Qb)-Vb{a^Q)]~l{k + k)Pa + ha''Pb 

+ -£s{aRa) + -V\aRab) ~ ^VaiaR) + uRa - I {3k - k) Sa'Sb - ka^Eb'Sc . (24) 

a a 2a 2 

The corresponding equations for Q and Qa are obtained from this by applying the Dirac duality transformations pOfl . 



III. EXACT SOLUTIONS ON A MINKOWSKI BACKGROUND 



In this section, we consider the Bianchi equations (O on the Minkowski background 
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where gAsdx'^dx^ = d-d^ + sin^ d(p'^ denotes the standard metric on 5^. A natural fohation is given by the shces 
t = const., for which Uadx'^ = —dt, although other foliations are possible. (In particular, it should be interesting 
to generalize the investigation below to hyperboloidal slices.) Since the spacetinie is spherically symmetric, it is 
natural to choose Sadx"- = dr. The corresponding vector field is defined everywhere except at the center r = 0. 
Furthermore, k = 2/r and kab = 0, so the evolution and constraint equations derived in the previous section simplify 
considerably. In this section, we assume that the source terms Rab and Sab vanish identically; however, we do not 
necessarily enforce the constraints Pa = Qa = since we are also interested in studying the propagation of constraint 
violations. 



A. Hcirmonic decomposition 



Since the spacetime is spherically symmetric and the equations are linear, it is convenient to expand the fields in 
spherical tensor harmonics. In the resulting equations, pieces belonging to different angular momentum numbers £ 
and TO decouple. Thus, it is sufficient to consider one fixed value of £ and m at a time. The decomposition of the 
fields Eab and Hab into spherical tensor harmonics reads 



E 


= -eo(i,r)F, 
r 




Ea 


= ei{t,r)VAY + fi{t,r)SA , 


Eab 


= 2re2{t,r) 


VaVb 


'^Y + 2rf2{t,r)V^ASB) 


H 


= -hoit,r)Y, 
r 




Ha 


= /ii(t,r)V,, 


Y + gi{t,r)SA , 


Hab 


= 2rh2{t,r) 


VaVb 


Y + 2rg2{t,r)V^ASB) 



(25) 



where Y = y^™(i?, (^) denotes the standard spherical harmonics, Sa = sa^^bY, and V denotes the covariant 
derivative on S"^. Similarly, the constraint variables can be written as 



P=-Poit,r)Y, PA = Pi{t,r)VAY + P2{t,r)SA , 
r 

Q = -Qoit, r)Y, Qa = Qiit, r)VAY + Q2{t, r)SA 



(26) 
(27) 



V2r 



The Newman-Penrose scalars are 

*o = 
*i = 

*2 = 
*3 = 
*4 = 

where rh^ = rm^. Using the identities 



-(e2 - g2)m^m''S7 A^ bY + - (/a + h2)m^m''V aSb 
r r 

^ {ei-gi)m^VAY ^ • u ^^A, 



V2r 



{h+hi)m''SA , 



— (eo - iho) Y, 
2r 



1 -A' 1 - A' 

(ei+5i)m \/aY y=-{fi - hi)rh Sa , 



V2r 



V2r 



-(e2 + g2)'m^m^^A^BY + -(/2 - h2)rh^rh^'^ aSb 
r r 



VaVi 



V^'V^aSb) 



*/ A - 

Y = --Way, 

-^Sa, 



(28) 
(29) 

(30) 
(31) 

(32) 

(33) 

(34) 



where A = (^— 1)(^+2), the Bianchi equations yield a set of two decoupled systems for the amplitudes (eo, ei, 62,51,52) 
(even parity sector) and {ho, /ii, /i2, /i, /2) (odd-parity sector). For £ > 2, the even parity sector is described by the 
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Pn — 


r 






CI 




A 




62 = 


-92 + 


1 

2r 




5i 




A 

2r 


3 


52 = 


-4 + 


1 





evolution system 

+ (35) 

(36) 
(37) 
(38) 
(39) 

where here and in the foUowing, a dot and a prime denote differentiation with respect to t and r, respectively. The 
evolution system (|35l3t)l37l38l39|l is subject to the constraints Po = Pi ^ Q2 = 0, where 

Fo = ^(^W-^^^ei, (40) 

Pi = ^(r'ei)'--e2--^eo, (41) 

Q2 = ^(^V)'~^52. (42) 

The odd-parity sector is obtained from this after the substitutions (60,61,62,51,52) ^ iho,hi,h2,—fi,—f2)- There- 
fore, it is sufficient to discuss the even parity sector, which is what we will do in the following. 



B. Exact solutions 



In this section, we discuss how to obtain exact analytic solutions to this constrained evolution system. We start 
with the special cases £ = and £ = 1 which, as we will see, are non-radiative. For £ — 0, the only equations are 
60 = and (r^eg)' = which yield the solution 

60 ~ — , M = const. (43) 

For £ = I, Eab and Hab vanish, and the evolution equations for 62 and 52, and the constraint equation Q2 = 0, 
are void. Taking a time derivative of the constraint Po = 0, and eliminating eo and ei using the evolution equations, 
one obtains 51 = c(t)/r^, where the function c is independent of r. The insertion of this information back into the 
evolution equations for 60 and 61 gives eo = — 26i -I- fc(r) for a function k which is independent of t. Substitution of 
these results into the evolution equation for 51 and the constraint Pi = gives c{t) = C2t + c\, k(r) = 2c2/r for some 
constants ci and C2, and 61 = C2/(2r) -I- f{t)/r^ for a function / which is independent of r. Finally, the insertion of all 
this into the evolution equation for eo yields / = C2t -|- ci. Thus, one obtains the most general solution in the sector 
£=l: 

C2 C2t^ + 2cit + 2co C2 ^ C2i^ -I- 2cit + 2co + Ci 

If we demand that the solution be stationary, then the corresponding odd-parity solution is 

ho = -^, '^i ^ ' /i = Oi J = const. (44) 

With the normalization y(^=o) = 1, — cosi?, Eqs. H43I44|) yield the linearized Weyl tensor belonging to the 

linearized Kerr metric 

-dt^ + dr^ + {dd^ + sin^ d dip^) + ^ {dt^ + dr^) - ^ sin^ d dipdt, 



where the linearization is performed about flat spacetime in the mass parameter M and the angular momentum 
parameter J. 
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Next, we consider the cases with £ > 2. We obtain two classes of solutions, one representing outgoing radiation 
and the other representing incoming radiation. If the constraints are satisfied, these solutions correspond to those 
presented in EDi IM| ■ To solve the above constrained evolution system, we first derive the constraint propagation 
system, which describes the propagation of constraint violations under the flux defined by the evolution equations 
(|6I7() . The constraint propagation system can be obtained cither by performing a multipolar decomposition of the 
evolution system H23I24|I with Rat — Sab — 0, or by taking a time derivative of Eqs. H40I41I42(I and using the evolution 
equations (|35l3(jl37l38l39|l . The result is 



Po 
Pi 



1 



Q2 = -7T P'l - -Pi 



Solutions to this system have the form 



Pn = — 



2r 



TT + rh'{r), 



Pi 



Hr) 



Q2 



where h(r) is a function of r only and where the function tt satisfies the master equation 

iii + iy 



2 



(45) 
(46) 
(47) 

(48) 
(49) 



with c = 1/2. Eq. (|49|l describes the evolution of constraint violations which propagate at half the speed of light. 
How to solve Eq. (|49|) will be explained below. 

Once the constraint variables Pq, Pi and Q2 have been obtained, we proceed as follows. First, using Eqs. (140141142(1 
we express ei, 62 and 32 in terms of Cq and gi and the constraint variables: 



l)ei = -cb'-rPo 



Ae2 = 



\g2 



-{r'9iy-rQ2 , 
r 



(50) 
(51) 
(52) 



where we have set 
equation for 0: 



^eo. Next, using these expressions in Eqs. H35|l and H38|l . we obtain the following wave 



- d't 



(53) 



where the source term S{t, r) depends on the constraint variables 7r(t, r) and h(r) and is given by 

£(1 + 1) 



S{t,r) 



[Srn' + 2n- 2rh{r)] ~ p/i'(r)]' 



Once Eq. (|53() has been solved for (j){t,r), the quantities ei, 62 and §2 are obtained from Eqs. H50I51I52|I . Therefore, 
the linearized equations reduce to the two master equations 



We now discuss how to obtain exact solutions to these equations. We start with the honiogeneous case where 
S{t, r) — 0. For the following, it is convenient to introduce for each i — 0,1, 2, ... the operators [50| 



ai = dr + - '^dr{r.), 
r 



and their formal adjoints 
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They satisfy the operator identities 



t t n2 ^(^+1) 

ae+ia'g^-^ = a'^at — —a^ H — 

As a consequence, for each ^ = 0, 1, 2, 

,2 , ^(^+1) 



(54) 



dt - di 



alal 



ala\_,...a\ [9^ _ 9^] 



(55) 



Therefore, in- and outgoing solutions to the hor nog eneous master equation can be constructed from in- and outgoing 
solutions of the one-dimensional wave equation [50|. For i = Q, the in- and outgoing solutions arc given, respectively, 
by (j)^fi{t, r) = Vb(r + t) and (f) y^.oit, r) = Uo{r — t), where Vq and Uq are smooth functions. For ^ > the solutions 
have the form 

(j)-\^^i{t,r) = a\a\_-^...a\VfXr + t), 
4)^,i{t,r) = a\a\_-^^...a\Ui>{r -t), 

where and Ug are sufficiently smooth functions. Explicit expressions for and (j) /^^i are given by 



K,£(t,r) = (-i)V 

0/,,(t,r) = (-l)V 



d 1 
dr r 

d_l 

dr r 



v,{r + t) = Y,{-iy 



^(2r)^-V«(r + t) 



3=0 



(^-J)!j! 
(2^-j)! 



(56) 
(57) 



where here and in the following, for a function F on the real line, i^^^-* denotes its j'th derivative. As an example, for 

i>^A^.r) = ^U2[r -t)~ -U!^\r - t) + u!^\r t). 

In- and outgoing solutions of the constraint propagation master equation (|49|l can be obtained in exactly the same 
way after replacing t by ct, ie., 



7r^/(t, r) = alal_j^...a\Wi{r + ct), 
■ny-^g{t,r) = a\a\_-^...a\Zi{r - ct), 



(58) 



for some sufficiently smooth functions Wg and Z^. 

Finally, we discuss the case S{t, r) ^ 0. Since we have already calculated the solutions to the homogeneous problem, 
it is sufficient to construct one particular solution of Eq. (|53ll . In the following, we assume that h=Q and that 7r(i, r) 
has the form TT{t,r) = a'gal_j^...a\Wg{t,r), with Wg a sufficiently smooth function of t and r. The latter assumption 
is no restriction of generality, since we can obtain W£ from tt by successive integration. 



Wi{t,r) = r / driri / dr2r2 



dri 



TT{t,ri) 



provided that 7r(i, .) falls off sufficiently rapidly as r ^ oo. To construct a particular solution (/'i(t, r) of Eq. (|53|l . we 
first notice that the operators 

Pin = 'irdr + m, 
m — ... — 2, —1, 0, 1, 2, ... satisfy the commutation relations 



Pma\ = a\pm-'i 
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As a consequence, 



If we make the ansatz 

4>i{t,r) = alal_j^...a\'ilj{t,r), 

and use relation ()55f) . we see that (pi is a particular solution if -0 satisfies the inhomogeneous one-dimensional wave 
equation 



With trivial initial data, this equation has the solution 

t r+t-T 



r-t+T 

Therefore, a particular solution of Eq. H53|) is given by 

t r+t-T 



<l>i{t,r) = ^■^-^^alal_^...a\ J J P2-uWi{t, s)dsdT. (59) 



r-t+T 



IV. SOLUTIONS OF THE INITIAL-BOUNDARY VALUE PROBLEM 



In this section, we analyze the Bianchi equations in the presence of artificial boundaries. Specifically, we solve the 
equations on a tubular subspace M = [0, T] x Br of Minkowski spacetime, where Bn is a ball of radius R in Euclidean 
space. Our goal is to impose boundary conditions on the timelike boundary T = [0, T] x dBp- which are perfectly 
absorbing in the following sense: for given initial data which is compactly supported in Bn and which represents a 
purely outgoing solution, the solution to the IBVP leads to the same solution as the solution to the global problem 
(without artificial boundaries). As discussed in the introduction, this turns out to be a challenging problem, even 
for simpler systems like the wave equation in more than one dimension 0, |3| . The strategy here will be to impose 
different boundary conditions on T which have been proposed in the literature, and construct exact solutions of the 
resulting IBVP by using the expressions derived in the previous section. With the help of these solutions, we analyze 
how "good" the boundary conditions are by looking at the amount of artificial refiection of constraint violating modes 
and gravitational radiation. This analysis enables us to construct different classes of boundary conditions, where each 
new class yields an improvement over the old one. 

We start in the next subsection with our crudest approximation for constructing outgoing boundary conditions, 
which consists of using the symmetric hyperbolic structure of the evolution equations (|6I7|I to freeze the incoming 
characteristic fields to their initial values. These boundary conditions yield a well posed IBVP. However, as we show, 
they introduce constraint-violating modes into the computational domain and therefore fail to be perfectly absorbing 
at a very fundamental level. 

In subsection IIV Bl first we specify CPBC which freeze the Weyl scalar '^q to its initial value. This is done either 
by using the method in [l^ . where suitable combinations of the constraints are added to the evolution equations, or 
by setting to zero the incoming constraint fields. Then, we show that the resulting boundary conditions introduce 
some reffections of (linearized) gravitational radiation. We quantify the amount of reflection by considering outgoing 
waves with wave number k and computing the reflection coefficient as a function of the dimcnsionless quantity kR. 

Finally, in subsection IIV CI we improve the boundary conditions considered in subsection IIV Bl In particular, for 
each L > 1, we give CPBC which are perfectly absorbing for outgoing gravitational radiation with angular momentum 
number £ < L. 



A. Freezing the incoming fields 



Let Sa be the unit outward normal one- form to the boundary dBn. For a symmetric hyperbolic evolution system 
of the form 
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with A^{u) = 1, the characteristic speeds and fields with respect to Sa are defined, respectively, as the eigenvalues 
and projections of u onto the corresponding eigenspaces of the matrix A'^{u)sa- For the evolution system H6I7|I . these 
are given by 

E = 2Re*2 , H = -2Im*2 , 

^Ea+ Sa'Hb = -V2i^iTna + §im,), 

V'ifc"'' = Eab + e(^a''Hi,)c = i^amafhb + ^orriamb), 











a 
2 






+ - 

^ 2 




A* = 


—a 




A* = 


+a 





where we use the same notation as in Sect. IIIBI and where /3° denotes the shift vector field. The ingoing fields are 
the ones with negative characteristic speeds (x. One way to obtain a well posed IBVP is to freeze the ingoing fields to 
their initial values !l^. In our choice of coordinates, with a — 1 and = 0, this boundary condition is equivalent to 
imposing 

dt^o ^ 0, dt^i = 0, 

where here and in the following we use the notation = to denote equalities which hold on the boundary dBn only. 

We analyze these boundary conditions using the harmonic decomposition of the previous section. As noted before, it 
is sufficient to consider the even parity sector. InviewofEqs. (|28I32() . we define the radial Weyl scalars -00 = 2(62—^2), 
■01 = ei — gi, 02 = eo, 03 = ei + gi, and 04 = 2(e2 + .92) which are functions of t and r only. Using the relations 
(|35l50l51l52fl . the master equation ((S2Jl, and Eq. (gSjl, we find 



BrTT + Srvr' + Gtt 
1)£(^+ l)(^ + 2)r4 ' i{e -!){£ + 2) 



^i = Z(7TTF+2' 

02 = ^ , (62) 

^3 = ^(7^-2' (63) 

6^0 — 8r7r + Srvr' + 67r 

" (£-l)£(^+l)(£ + 2)r4 + 4(^ -!)(£ + 2) ' ^ ^ 

where we have introduced the operators b± = r^{dt T dr) and again assumed that h = Q for simplicity. In what 
follows, the expressions 

{b^r^^At.r) = 1^ 1: (-1)^+- f 1?]),^., (20^--^^/^"^(^ + t\ m^0,l,2,...£, (65) 

(6_)™0^,(.,.) . (-1)V""E ' (2^+"^^ ■)! i)!.! ^'^^'"'"''^"'^^ + '^' m . ^ + 1, £ + 2, (66) 

(6_)™0/,,(t,r) = ^(-l)'^" r. ■, {'^ry-'ul'\r-t), m = 0,l,2...£, (67) 

(6_)"0/,£(i,r) = 0, m-^+1,^ + 2,... (68) 

will be useful. They can be derived from the explicit expressions (|56I57|) by induction in m. Corresponding expressions 
for (&+)™0 can be obtained by fiipping the sign of t and interchanging 0x,^ 4' We see from these expressions 
that if TT = and = then along the outgoing null rays r — t = const., we have 6_0 = 0{r'^), b'^_4> = 0{r^) and 
6+0 = 0{P), b'^cf) = 0{r^). Therefore, the radial Weyl scalars obey 

iP^ = Oir"^'^), r-t = const. (69) 
for s = 0, 1, 2, 3, 4. This is consistent with the peeling theorem [s^ . 



14 



Next, we construct exact solutions of the IB VP described by the evolution equations (|49I53() . initial data for tt and 
(j) on the interval (0, R), and the boundary conditions dt'^^o = dt'^fi = imposed on a sphere of radius R > 0. These 
solutions have the property of satisfying the constraints initially {ie. 7r(0,r) = and Tr{0,r) = for all < r < i?), 
but violating the constraints at later times. For the sake of avoiding unnecessary complications, we restrict ourselves 
to the case £ — 2, although one should be able to construct similar solutions for higher i. 

To construct these solutions, we start with a smooth function F : (0, oo) R which is zero on the interval (0, R) 
but non-zero for r > R, and set 

7r(t,r) = alalF^'^^r + ct), 

where F^^^ denotes the fifth derivative of F. By construction, this solves the constraint master equation (|49|l . and 
7r(0, r) = 0, Tr(0, r) = for all < r < R. Choosing h = guarantees that the constraint variables Pq, Pi and Q2 
have trivial initial data as well. Next, using the construction procedure outlined in the previous section, we obtain 
the general solution of the inhomogeneous wave equation (|53|l . The result is 



0(t,r) 



(l)y{r-t)+(j)'K^{r + t) + 2AF'-^\r + ct) - 6rF'-'^\r + ct) 



(70) 



where y and are (up to this point) arbitrary smooth functions. In order to determine these functions, we insert 
the general solution 17U|I into the boundary conditions dt'^Q = dt'i'i = 0. Using the expressions H6(JI61|I . we obtain 



4 2 

4rVo (t, r) = cl>y{r-t) + 4>^{r + t)- 2r(f>i^ {r + t) + 2r^(l)^^ [r + t)-- r^(j)^^ [r + t) + - r^^^ {r + t) 



6P(^) (r + ct) - QrF'^^^ (r + ct) + Gr^P^^) _^ ^i) - 2r^F'^^^ (r + ct) 



-rVi(i,r) = (t)y{r - t) ~ ^r(l>y{r - t) + ^f-^ir + t) - ^ r^l^^ (r + t) + r^^^^^ (r + i) - |r^0^'(r + i) 
+ 24P(3) (r + ct) - aOrP^'^) (r + ct) + y r^P^^) (r + ct) - 5r^F^^^ {r + ct) + r'^F^'^^ {r + ct). 
The combination B = r'^tpi + Ar'^tpQ + 2r^tl'o gives 

B{t,r) = i r^cjj^^ir + t) - ^ r^F'-^^r + ct) + 3r^F<^^\r + ct) ~ 3r^F^'^\r + ct). 



.2^(2)/ 



.3^,(3), 



(71) 



(72) 



(73) 



Therefore, the boundary conditions ?/'o(t, P) = ?/'i(t, P) — for all t > imply B{t,R) — for all t > 0. After 
integrating Eq. H73|) and setting five integration constants to zero, this condition yields 

cf)^(R + t) = ^ P(P + ct) - 2PP(i) (P + ct) + 2p2p(2) (i? + ct) 

for all t > R, thus determining on the interval (P,cxd). On (0, P] wc simply set (f>'\ to zero which means that 
initially, the solution does not contain any ingoing radiation. Plugging this and the expression H71|l into the boundary 
condition "00 (t, P) = for alH > fixes 4)y{R — t) for all t > 0. The final result is 



4>y{r-t) = 



P3 



-72P(z) + 216PP(^Hz) - 324p2p(2)(z) + 216P^p(^)(z) 
- 81P'^p(^)(z) + 18P^p(^)(z) - 2P^p(6)(z) 



z=c(3-R-r+t) 



t) = — ^ [p(z)-2PP(i)(z) + 2p2p(2)(z)' 



c{R+r+t) 



(74) 
(75) 



Therefore, we have constructed explicit solutions which satisfy the boundary conditions dt'^o = dt'^i = obtained 
by freezing the incoming characteristic fields of the symmetric hyperbolic system H6I7() . These solutions have the 
property that the constraints Pa — Qa — are satisfied initially, but violated for t > 0, thus providing an explicit 
example which shows that freezing the incoming characteristic fields to their initial values at the boundary is not 
always compatible with constraint propagation. This fact has also been observed in numerical simulations |l6lll7| . 



B. Constraint-preserving boundary conditions: Freezing ^0 



Here, we improve the boundary conditions considered in the previous section. Our goal is to formulate the evolution 
problem in such a way that solutions belonging to constraint-satisfying initial data automatically satisfy the constraints 
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everywhere on Bfi and at each time t > 0. There are two ways to achieve this. The first approach T^l modifies the 
evolution equations by adding suitable combinations of the constraint equations to them in such a way that the 
resulting constraint propagation system is symmetric hyperbolic and does not contain any normal derivatives at the 
boundary. Consequently, the constraint-preserving property of the boundary conditions is automatic. The second 
approach leaves the evolution equations unchanged, but replaces the boundary condition = with a carefully 

chosen boundary condition which guarantees constraint propagation. 
In the first approach, one chooses 

Rab = +S(a£b)''''scQd , (76) 
Sab = -S(„eb)=^S,Pd , (77) 

instead of Rab = Sab = in the evolution equations (|6I7|I . where Pa and Qa are given by Eqs. © and 0, respectively. 
One can verify that the resulting evolution system is symmetrizable hyperbolic and that the characteristic speeds and 
fields with respect to Sa are unchanged, with the exception of one important difference. The presence of the term 
£sEa in Pa (see Eq. (|20f) ') cancels the corresponding term in the evolution equation for Ha- Similarly, the term £sHa 
is canceled in the evolution equation for Ea- This changes the speeds of the characteristic fields V"i~'' and V^"*^' from 
±a/2 — l3°'Sa to —(3°'Sa. For our coordinate choice, with = 0, 5*1 is no longer an incoming field, and therefore, 
no longer requires boundary data. The only remaining boundary condition is the one involving '^q. The definitions 
(I76|l and (|77|l yield a different constraint propagation system than discussed in the previous subsection. Using Eqs. 
(I23I24|I . we obtain 

£nP = e"'^- ("'Qb) + ^ " 3fc) P, (78) 



£nPa = -T^Ea'Vr^ia^Q) - kabS^'Qc + ( — - ) Sa^Qb - kPa + 2ka' Pb • (79) 

The corresponding equations for Q and Qa are obtained from this by applying the Dirac duality transformations H10|l . 
For the case of linearization about Minkowski space in the natural foliation where a = 1, Z?" = 0, kah = 0, k = 2/r, 
and kab — 0, we obtain, using the harmonic decomposition H2t)l27|) . 

Po = -^^^^Q2 , (80) 
r 

Pi - -Q2, (81) 
r 

Q2 = -^(2Pi+Po)- (82) 
2r 

These equations are ordinary differential equations in time. Therefore, initial data which satisfies the constraints 
automatically yield constraint-satisfying solutions. 

In the second approach, one sets Rab = Sab = as before, but replaces the freezing boundary condition dt'^i = 
with CPBC which can be constructed as follows. For Rat = Sab = 0, the constraint propagation system H23I24|I is a 
symmetric hyperbolic system whose characteristic speeds and fields with respect to the radial field Sa are 

M = -r-s, , P, Q, (83) 

f^^-!l-f3-Sa, Wt^ =Pa+Sa''Qb, (84) 
P^Sa , Wi+^ ^Pa~ Ea^'Qb ■ (85) 

If = (or more generally, if l3°'Sa < 0), the homogeneous boundary condition 

W^-^ = 

guarantees that the unique solution to the constraint propagation system with zero initial data is zero. In terms of 
the constraint variables h{r) and 7r(f,r) introduced in Eq. 1481) . this boundary condition yields 

-tt + tt' ^ -h{R), 
c c 

guaranteeing that solutions of the constraint master equation (|49|l which satisfy h{R) = 0, 7r(0, r) — 0, and 7t(0, r) — 
for < r < R are trivial. Therefore, there cannot exist solutions with constraint-satisfying initial data which violate 
the constraints at some t > 0, as we encountered in the previous subsection. 
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From now on, we assume that the constraints Pq = Pi = Q2 = are exactly satisfied, and analyze solutions to the 
Bianchi equations on [0, T] x Bf/ which satisfy the boundary condition 



dtiJo = 0. 



(86) 



As indicated above, ipo does not vanish exactly at a finite radius for the purely outgoing solutions y^g, but falls off 
as r"* on the outgoing null rays r — t = const. Therefore, imposing the boundary condition l|86() at finite radius 
r ^ R < 00 yields reflections of gravitational radiation. In other words, solutions to the IBVP will consist of a 
superposition of a purely outgoing and a purely ingoing solution. In order to quantify the amount of reflection, we 
first consider monochromatic quadrupolar waves of the form 



(87) 



where fc is a given wave number which is assumed to be different from zero and 7 a (yet unknown) amplitude reflection 
coefficient. Introducing this ansatz into the boundary condition (|86|l yields 



1+7 



l + 2ikR-2{kRf 



^^{kRf + l^kRf 



As is easy to verify, the expression inside the bracket is never zero. Therefore, we can solve this equation for 7. 
amount of reflection is given by 



l72(fci?)| 



l-^{kRf + -{kRf 



-1/2 



The 



(89) 



where the subindex 2 refers to the fact that we are considering quadrupolar waves. The reflection coefficient |72(fci?)| 
versus kR/2 is plotted in Figure^ There is a global maximum at kR — •y/3/2 where |72(fci?)| — 2. For kR ^ 1, 
|72(fci?)| decays as (kR)^^. Therefore, the boundary conditions are very accurate provided the size of the domain 
is much larger than the characteristic wavelength of the problem. On the other hand, if the size of the domain is 
comparable to the characteristic wavelength, the reflection coefficient is of the order of unity. How to improve this 
boundary condition is explained in the next subsection. 

Using the expressions (|65I67() . the above analysis can be repeated for arbitrary £ >2. The result is 



\MkR)\ 



where the polynomials pf^rn{z), \m\ < £, are given by 



pe-2{-ikR) 



Pi,2{ikR) 



(90) 



Pl.Tuiz) 



£+m 

E 

3=0 



m)!(2£- j)! 



-J)!j! 



i2zy 



(91) 



The refiection coefficients qi — |7^(fci?)| versus kR/£ for different values of £ are shown in Figure ^ It can be seen 
that qi is of the order of unity for kR/£ < 1 while for kR/£ » 1, decays very rapidly. From Eq. it follows 
that for large kR, |7^(fci?)| decays as {kR)~'^. Although for fixed kR, the reflection coefficient gets larger when £ is 
increased, this is not an issue for most physically interesting scenarios, since the first few multipoles usually dominate. 
In particular, if the solution is smooth, amplitudes corresponding to different values of £ decay rapidly as £ 00. 
Therefore, even though for high £'s the reflection coefficient is large, it does not introduce a large overall error since 
the corresponding amplitudes of the solutions should be very small. 

Figure |21 shows in more detail the amount of reflection if the outer boundary is placed at a few multiples of the 
characteristic wavelength of the problem. Clearly, this amount of reflection is very small (0.1% or less for R greater 
than or equal to one wavelength and £ — 2, and less than 0.0065% for R greater than or equal to two wavelengths 
and £ = 2). 
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kR/l 

FIG. 1: Reflection coefficient qi = \'yt{kR)\ as a function of kR/d (I = 2, 3, 4, 5, 20), for the boundary condition dt'^o = 0. 











1=2 

1=3 

1=4 

1=5 

1=20 












R = 1 wavel 


ength, 1 = 3 












\ \L 
\ w 








\ V 

\ ^ 








\ >L 


R = 1 wa\ 


felength, 1 = 2 

















— 7e-05 



4 

kR/l 




FIG. 2: Close-up of regions 2 < kR/£ < 6 and 6 < kR/£ < 10. 



C. Improved constraint-preserving boundary conditions 



As we have seen in the previous subsection, the boundary condition H86|l is not perfectly absorbing. If the outer 
boundary is a sphere of radius R, and for monochromatic waves with wave number fc, there are reflections where the 
reflection coefflcient is proportional to (kR)^'^ for large kR. Although these reflections can be made arbitrarily small if 
the boundary is pushed sufficiently far away, there is significant motivation for improving the boundary conditions. In 
particular, it may not always be possible to push the outer boundary far into the wave zone in numerical simulations, 
especially for those in three space dimensions, because of the high computational cost. Even if this can be achieved, 
it may still be desirable to decrease the artificial reflection in order to achieve better accuracy. 

Our goal here is to find boundary conditions which are perfectly absorbing at least for all multipoles i = 2,3, ...L, 
where L is a given maximum. This means that for initial data which is supported on the interval (0, R) and which 
corresponds to a purely outgoing solution (j) y^^iit, r), the solution of the IB VP for i > is uniquely given by (f) /'.^(t, r). 

One way to achieve this goal is to rely on the identities H54|l and the fact that ^/'z solves the homogeneous master 
equation H53() for each £ gN. Using this approach, we find 

aia2...ae(j>/',e = aia2...aial(j)y^^e-i = -aia2...ai^idt<j)^^e-i = ■■■ = {-lY{dtf'^<j>^fi{r - t). 



This expression vanishes identically if we apply the operator 6_ = r^{dt 
for our perfectly absorbing boundary condition on the field (f> = r'^eo is 

b-aia2...ae(t> = 0. 



dr) to both sides. Therefore, a candidate 



(92) 
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However, a problem with this condition is that it is only quasi-local in the sense that it is different for each £. Therefore, 
a numerical implementation of the IBVP requires performing a harmonic decomposition of the electric and magnetic 
fields Eab and Hab near the outer boundary so that (f) can be computed and the boundary condition applied. 
An alternative approach is based on the observation that for all £ > 2, the outgoing solutions (p y^g satisfy 

(6_)^+V^,£(i,r) = 0, (93) 

which follows directly from Eq. (|68|l . We therefore impose the boundary condition (b^)^^^(f) ^ 0. It turns out that 
this boundary condition agrees precisely with the hierarchy of conditions given in for the flat wave equation in 
three space dimensions. There, it was also shown that the boundary conditions yield a well posed IBVP for the wave 
equation and that the error with respect to the solution on the unbounded domain (measured in an appropriate norm) 
decays as as the radius R of the outer boundary goes to infinity. In order to allow for a static contribution 

to (j>, we impose the boundary condition 

(5_)^+iat0 = O. (94) 

In the appendix, we prove by deriving a suitable estimate that the resulting IBVP is stable, and that the initial 
data uniquely determine the solutions. As a consequence, the boundary condition H94|l is perfectly absorbing for all 
multipolar waves with £ < L. In view of Eq. H60() this boundary condition is equivalent to the condition 

(6_)^-i(r4atV^o)^0 (95) 

on the radial Weyl scalar ipo, provided that L > 1 and that the constraints are satisfied. Therefore, for L > 1, the 
boundary conditions (|94|l can be reformulated as boundary conditions on the incoming characteristic field ^'o and its 
derivatives. This sheds some light onto the meaning of the freezing boundary condition: it is the first member 
of a sequence of boundary conditions with increasing order of accuracy. By construction, the boundary condition 
H95|l is exactly satisfied for all outgoing linear gravitational waves with £ < L. The uniqueness result in the appendix 
also implies that it sets to zero any incoming gravitational radiation. Furthermore, the boundary conditions l|95|) 
is local in the sense that it does not depend on £. Thus, a numerical implementation does not require a multipolar 
decomposition. 

Finally, we compute the amount of artificial reflections for solutions with £ > L. In order to do so, we generalize 
the ansatz Eq. (|57|l to arbitrary £: 

I/, \ t t / ik(r-t) I -ik(r+t)\ 

0(r, r) = a^...a| I e ^ ^ + 7e )• 
Inserting this into Eq. 195|l . using Eqs. (|6()I65I67(I . and assuming that the constraint variable tt is zero, we obtain 

l7L,£(fci?)| = : ^>L, (96) 

where the polynomials pi,m{z) are given in Eq. (jni. In particular, |7L,f(fci?)| falls off as (fci?)-2(i+i) for large kR. 



V. EFFECTS DUE TO THE BACKSCATTERING 



In this section, we want to analyze how the results obtained in Sects. IIVBI and IIV CI are modified if instead of 
considering linear wave propagation on a flat spacetime, the background is curved. In particular, we are interested in 
a physical situation involving a localized region of space where strong gravitational interactions take place, and where 
outside this region, the gravitational field decays rapidly to flat space. Therefore, far from the strong field region, 
we can expect spacetime to be accurately described by a perturbed Schwarzschild metric of mass M representing 
the total mass of the system. We place a spherical boundary of radius r — R, where r is the area radius of the 
Schwarzschild background, and assume that 2M/R <C 1. In the following, we generalize the constructed in- and 
outgoing wave solutions to include first order corrections in 2M/R, and then compute the first order correction terms 
to the reflection coefficients found in the previous section. For simplicity, and since we are only interested in the 
qualitative behavior of the correction terms, we restrict ourselves to perturbations with odd parity. The effects of 
second order corrections in 2M / R, corrections due to J / R^ , where J is the total angular momentum of the system, 
and corrections emanating from nonlinearities (see Ref. ^53J for an estimate on the errors introduced by neglecting 
the nonlinearities of the theory) will be considered elsewhere. 
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A. Odd-parity linear fluctuations and derivation of a master equation for Im5'2 



As shown in Ref. |54j , linear odd-parity metric perturbations about a Schwarzschild black hole can be described by 
a master equation for Im54'2, where 6'i>2 denotes the linearization of ^'2- Since lm^'2 is a scalar field that vanishes on 
a spherically symmetric background with an adapted Newman-Penrose null tetrad, its p erturbation is invariant with 
respect to infinitesimal coordinate transformations. Additionally, one can also show |5j| that lm5^2 is invariant with 
respect to infinitesimal rotations of the null tetrad, and is therefore well-suited for describing odd-parity grav itational 
perturbations. It turns out that the master equation for this quantity is the Regge- Wheeler equation |5fi| . In this 
subsection, we briefly review the derivation of the Regge- Wheeler equation for Im^'2. For simplicity, we assume that 
the background is written in standard Schwarzschild coordinates (t, r, z?, ip) for which 

n^da = -dt , s^da = adr , labdx^'dx^ = idd^ + sin^ ^ dip^) , 
a 

_ - _ n - 

r 



where a = ^^1 — 2M/r. The corresponding electric and magnetic parts of the Weyl tensor are 

M 

Eab ^ -^{jab - 2SaSb), Hab = 0, 

where the circles on the top of Eab and Hab indicate that they are background quantities. Linearizing the evolution 
and constraint equations (|6I7I8I9|I about this background yields the system 

(97) 
(98) 
(99) 
(100) 

where n, Saba D, and a'^ refer to the background geometry, and where Eab and Hab denote the parts of the perturbed 
electric and magnetic components of the Weyl tensor which are trace-free with respect to the background metric 
hab = SaSb + Jab- Also, indiccs are raised and lowered with the background metric. The source terms Rab, Sab, Pa, 
and Qa depend on the perturbations of the shift, SP"", the perturbations of the metric, dhab, and its first spatial 
derivatives, and the perturbations of the extrinsic curvature, Skab- Performing a change of infinitesimal coordinates 
if necessary, we can obtain (5/3° — 0. In this case, we find that the source terms Rab, Sab, Pa, and Qa are given by 

1. -rf.. \ _^ 









= Rab 




- e^diaiD-' - 


f 2a')E\^ 


— Sab 






D^'Eab 


= Pa , 






D^Hab 


— Qa 



Rab = 5[ E^a"5kb)c - ^habE Sk^d ) - 2Eabh^''5kcd , (101) 



C ^ cd 

iJab — ~fc(a 



Set \ 

Shij^ca'^Ede + C'^b)cEde + 2E},)cDd ( — 



(102) 



Pa = {D''Ea-^DaE"yhb, + Eabh^''C\d+\Ek'C\a, (103) 
Qa = -ea'"'Eb''Skcd , (104) 

where 

C\b = \h"' {DaShtd + DbShad - DJhab) 

are the linearized Christoffel symbols. Here, we have also used the background equations 

£nShab = 2Skab , e,d{a{D^ + 2a^)E\) = 0, D^'Eab = 0, 



which imply 



£cda{P>'^ + 2a'') E b — EabcO-dE ■ 



Next, we perform a 2 -|- 1 split of Eqs. (|97I98I99I100|I as described in Sect. UTOl Notice that the 2 + 1 split is 
with respect to the unperturbed Schwarzschild metric, for which the assumptions made below Eq. H15|l on the vector 
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fields 3°" and n° hold, and not with respect to the perturbed Schwarzschild metric. Using the odd-parity sector of 
the harmonic decomposition 125|) for Eab and Hat, and a corresponding odd-parity decomposition for 5hab and Skat, 
namely, 

6hab 2<T(t,r)s(aS'b) + 2ri^(t, r)V(aS'b) , 
Skab = 2Trait, r)s(^aSb) +2rTr^(t,r)\'(^aSb) , 
we obtain the following equations for ^ > 2: 



ho £{£+!) 



a 2a 



a r 
2r ~" ' 4r' 



fi = - 



£{£+!) M 



5M 



M - 



hi = 0, 



a 2 7 \/ 7 1 7 3Af 
^(r /ii) 112-— hQ = 

^(rVi) - -/2 





7r<j , 






M 






A 






(7) 


V 








r 



(105) 
(106) 
(107) 
(108) 
(109) 



A master equation for (f) = r^/io is obtained as follows. First, we use Eqs. (|108|l and p07|l in order to eliminate 
/i2 and hi in Eq. H106(l . Then, we use Eq. H105I) and the definition of the extrinsic curvature, a — 2aiTa., in order to 
eliminate /i from the resulting equation. This yields 



-DM TTo 



(110) 



a r" 

To get an equation for <j) alone, we need a relation between and (j). This is obtained by linearizing the equation 

Hab = —£cd{aD'^k'^b), 

which expresses the magnetic part of the Weyl tensor in terms of the curl of the extrinsic curvature. This yields 



^ = -^(^+l)r7r<, , 

and leads to the Regge- Wheeler equation [s^ 



0(t,r) =0 



(111) 



(112) 



for (j). As in the flat spacetime case, the linearized Newman-Penrose scalar is entirely determined by (j). To see 
this, we first linearize Eq. H11() '8T| and obtain 



-00 



m^m^VA^B , Vo = 2 ( /i2 + /: 



M 



Using Eqs. (|105I109|I and Eq. pil|) . we re-express the above expression in terms of alone, giving 



(£- 1)^(£-H l)(^ + 2)r4 ' 



(113) 



where &_ = r^(a "^dt -\- dr). This generalizes Eq. (|60|l to a Schwarzschild background. 



B. Construction of in- and outgoing wave solutions to first order in M/R 

Next, we generalize the in- and outgoing wave solutions constructed in Sect. IIII BI to the case M ^ 0. This means 
that we have to solve the new master equation ()112|l . Since we are only interested in cases where M/r <C 1, it is 
reasonable to expand the equation in factors of M/r and to consider the first order corrections in M/r only. This 
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expansion might depend on the chosen coordinates. For the following, it is convenient to introduce the tortoise 
coordinate which is defined by 



r*(r) = 



ds 



1 - 



2M 



= r-AM + 2Mlog (— 1 



AM 



Using this, we can rewrite the Regge- Wheeler equation as 



- dl 



2M 

r 



f£{i+l) 6M\ 



cf){t,r) = 0. 



(114) 



Therefore, if r is very large compared to I and M, in- and outgoing solutions are, approximately, given by (t>\^{t, r) w 
+t) and (j) « U{r^,—t), respectively. Notice that is not analytic in 2M/r at 2M/r = 0, so it is not clear 

if, for example, (l)/'{t,r) ^ U{r — t) (r* replaced by r) is a good approximation for the behavior of outgoing solutions 
in the asymptotic regime. For this reason, it seems more appropriate to use the coordinates (t,r*) to describe the 
asymptotic behavior of the solutions. On the other hand, the potential term appearing in Eq. (|114ll is not analytic in 
2M/r* at 2M/r* = 0, so we cannot expand it in terms of powers of 2M/r* near 2Af/r* — 0. In order to circumvent 
this problem, we introduce the new coordinates 

T = t + r — r.^, , P — f, 

in which the Regge- Wheeler equation can be written as 



di~dt 



where the operator B is defined by 



dr + dp 



T, P), 



(115) 



If we neglect the right-hand side, this equation reduces to the fiat space master equation which has the outgoing 
solutions (f) y^iij^p) constructed in Sect. IIII Bl These outgoing solutions have the correct asymptotic behavior since 

(j) y/ir, p) « uf\p — t) = u\^'^ (r* — t). We also see that for these solutions, Bcj) y^i{T, p) decays as p^^ — r^^, so the 
right-hand side of (|115ll is small. Therefore, given R > 2Af , we expect that we can write the solution in terms of an 
expansion in powers of 2M / R as 



{t,p) = aeipY ai-i{p)'' ...aiip^Uip - t) ( ) 3k{T,p), 



(116) 



k=l 



for all /9 in a neighborhood of R, where here and in the following, aeXpY = —dp + i/ p. In Ref. |43j |. a similar expansion 
was used to obtain solutions of the Teukolsky equation 44] on a Schwarzschild background, and was shown to converge 
absolutely. Plugging the expansion (|116|) into Eq. H115|) yields the following hierarchy of partial differential equations: 



£{£+!) 



R 

9k{T,p) = Bgk-i{T,p), 

P 



k = 1,2,3, 



(117) 



where go{T, p) = ai{py ...ai{pyU{p ~ r). In Sect. IIII Bl we learned how to solve such equations using integral 
representations of the solution operator of the fiat wave equation. 

In the following, we give the explicit solution for the first order correction (fc = 1) of quadrupolar waves (£ — 2). 
The solution can be written as 



51 (t, p) 



where the integral kernel K2 is given by 



3R 



L/W(p- 



K2{t,p,x) = a2{p)^ai{py 



R 
1 



K2{t, p,x)U{x)dx, 



(118) 



p-T 



[T + p + xy 



3 r - 



2w' 



2u>' 



2p 



X > p 
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and satisfies 



K2{t,p,x) = 0, 



K2{t,p,p-t) = ^ 



15 



(119) 



{dr +dp)K2iT,p,p-T) = 



30 

75 



Notice that for t > and p > 0, the function x i— s- ii'2(''', Pi x) is bounded from above by the function 

30 1 



X ^ Ml (x) 



p2 (p + xy 



on the open interval x > p — t. Therefore, if U is continuous, supported on the interval (0,00), and bounded, then 
the integral in Eq. pi8|l exists for all r > and all p > 0. Using the properties H119|l . it is not difficult to verify 
that gi indeed solves Eq. H117|) for ^ = 2 and fc = 1. Notice that if U is supported in [ri,r2], where < ri < r2, 
the zeroth order solution go{T,p) = U^'^\p — t) — 3C/(^^(p — r)/p + 3C/(p — t)/ p^ is supported in [ri + T,r2 + r] for 
each r > 0. In particular, for each fixed pi > 0, go{TiPi) vanishes for r large enough. This is a manifestation of 
Huygens' principle which holds for the fiat wave equation in odd space dimensions. The first order correction term 
gi vanishes for p > r2 + r, but not necessarily for p < ri + r. This is the effect of the backscattering. Nevertheless, 
for each fixed pi > 0, o i (r, pi ) converges to zero as r — > cx3. This can be shown by using Lebesgue's dominated 
convergence theorem|82l| and noticing that the function x 1— > i^2(''', pi, a;)C^(a;), which is bounded by the integrable 
function x i— > Mi{x)\U{x)\ on the interval x > pi — r, converges pointwise to zero as r ^ 00. 

Summarizing, we have obtained outgoing, approximate solutions of the Regge- Wheeler equation for 1 — 2: 



4>^{t,r) = [/(2)(r,-t)--C/(i)(r,-t) + 4{/( 



2M 



4r2 ^ ^4 



- t) 

K2{t + r — r^,r,x)U{x)dx 



O 



/2M 



Since the Regge- Wheeler equation is time-symmetric, corresponding ingoing solutions are obtained from this by merely 
fiipping the sign of t: 



2M 



K2{—t + r — r^,,r, x)V{x)dx 



O 



f2M 



Using Eq. (|113|) and the fact that 6_ = a ^r^{dt + 9^, ), we compute the corresponding expressions for -0o: 



U{r^ -t) + 



2M 



'2U{n ~t) + -C/^^Hn -t) + - I k{l + y)U{r, - t + 2ry)dy 



O 



(120) 



^ox(i,r-) = (vir, + t) - 2rV<-'^{r, +t) + 2r^V<-^\r, + t) - ^r^V^^^r, + t) + Ir^V^^Hr, + t) 

4q,^j,4 \^ 3 3 



2M 

r 



00 

1 2t.(2). , ,N 1 3t.(3)/ , .X , 1 f V^r^ + t + 2ijj)dy 



O 



f2M 



(121) 



where k{w) = hw^^ + 4w"^ -I- 3w~** + 2^^^ + w^"^ . Taking into account the fact that ot"^ = 1 + 2Af/r + 0(2M/r)^ 
and replacing U{x) by G{x) = U{~x)^MU'-^\-x)l2, the result for a-'^^^^/^^ agrees with Eq. (4.18) of Ref. [ig. 



23 



C. Reflection coefficient for the boundary condition dt'^o = 



In order to quantify the amount of artificial reflections caused by a spherical artificial outer boundary at i? S> 2M, 
we consider as before monochromatic waves of the form 



f/(r, -t) = V{r, +t)= je-^'^'-'-'+'K 



(122) 



where 7 is the amplitude reflection coefficient. Introducing these expressions into Eqs. H12UI121|I and setting 
dttpo.y{t^ R) + dti'o,\{t, R) — 0, we obtain the result 



72 kR. 



2M\ 

~r) 



l72(fci?)| 



l + —E{kR) + 0[—- 



(123) 



where |72(fci?)| is the reflection coefficient given in H89|l (which is valid for M = 0), and the function E{z) is given by 

1 



i?(z) = --z6|72(z)P 



8z^ - 13 - {2z^ - 4) / fc(l + y) cos{2zy)dy 



(124) 



In deriving this result, we have used the integrals 

oc 

cos(2zy) 



Cniz) = 



(1 + yY 



dy, Sniz) 



sin(2zy) 
(1 + 2/)" 



dy, 



and the relations 



for n > 2, which imply that 



1 2z 

C„+i(z) = - [1 - 2zSn{z)] , Sn+l{z) = C„(z), 

n n 



lim {2zfCn{z) = n, lim (2z)S'„(z) = 1, 



for all n > 2 and 

00 

4 - z2 + - 2z2 + C2{z) + Q ^3 - 2zj S2{z) = y + y) cos{2zy)dy. 



Since E{z) —2 as 2; ^ 00 it follows that \j2{kR,2M/ R)\ still decays as {kR)~^ for large kR. In fact, for kR 
sufiiciently large, the reflection coefficient is smaller than the corresponding flat space coefficient provided that 2M / R 
is small enough. 



D. Reflection coeflScient for the improved boundary condition 



Finally, we repeat the above analysis for the boundary condition 

{dt+dr){r^dti^o){t,R) = {), 
which is perfectly absorbing for M — (see Eq. (|95f) '). From Eqs. (|120I121|I we first obtain 



2M 



C/(n -t) +r[/(i'(r, -t) - 15 



U{r^ ~ t + 2ry)dy 



4r(a* + 9,)(rVov)(i,r) = K^'V^^) [r, + t) + O {^-^ 



(125) 



O 



/2Af 
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2M/R 0.1 




FIG. 3: Reflection coeflicient 52 = \"^2(kR,2M / R)\ truncated to first order in 2M/R as a function of kR and 2M/R, for the 
boundary condition dt'^o = 0. 



2M/R 




FIG. 4: Reflection coefficient q2 = \'^2{kR,2M / R)\ truncated to first order in 2M/R as a function of kR and 2M/R, for the 
boundary condition dt'^Q = 0. Showing surface for 3 < kR < 8. 



Using the monochromatic ansatz H122fl . we obtain 

2M 



72,2 kR, 



R 



2M~,, ^ /2M 



(126) 
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where 



E{z) = Rl - lhCj{z)f + {z~ lbSj{z)f 
Az^ I 



1/2 



(127) 



For kR ^ 1, the reflection coefBcicnt goes as (2A//i?)(fci?)^^. Because of the presence of the smah factor {2M/R), 
there is a significant improvement over the boundary condition dt'^a = 0. In Figurel^l we plot the ratio E{kR) /\j2ikR)\ 
as a function of kR. This plot, together with the asymptotic expansion 2E{z) /\j2{z)\ = 1 — Sz^^ + 0{z~^), suggest 
that for kR > 1.04, this ratio does not exceed 0.5. Thus, we conclude that with corrections for backscatter, our 
improved boundary condition gives a reflection coefficient which is M/R times smaller than the one for the freezing 
condition for kR > 1.04. 




FIG. 5: E{kR)/\-i2{kR)\ versus kR. 



VI. CONCLUSIONS 

Numerical relativity gro ups around the world have begun to calculate binary black hole merger waveforms 

with the goal of providing waveform templates for the detection and 
interpretation of gravitational wave signals from instruments such as LIGOfs^, VIRGO'sB'l and LISA'S^. To be 
useful templates, the calculated waveforms must be as accurate as possible. In particular, if numerical binary black 
hole simulations are performed on a finite computational grid with an artificial outer boundary, it is critical that this 
boundary be as seamless an interface as possible between the physical scenario and the computational grid. Towards 
this end, we have constructed a hierarchy Sl, L = 1, 2, 3, ... of boundary conditions which are perfectly absorbing for 
linearized waves with arbitrary angular momentum number € < L on a Minkowski background with a spherical outer 
boundary. For a nonlinear Cauchy formulation of Einstein's vacuum field equations, these boundary conditions can be 
formulated as follows. Let t be the time-like coordinate compatible with the foliation by space-like hypersurfaces 
(ie. such that Tit = {t = const}), and let r be a radial coordinate which has the property that the two-surfaces St^r 
of constant t and r are approximate metric spheres with area Aur"^ for large r. We assume that the outer boundary is 
described, for each t > 0, by the two-surface 5'*,^. Let n° be the future-directed unit normal to the surfaces and 
let s°' be the normal to the surfaces St.r tangent to Et. Finally, let v"" and be two mutually orthogonal unit vector 
fields which are normal to n°- and s°, and define the real null vector Z° = (n" -I- s°')/\/2 and the complex null vector 
to" = (w" -I- iw°')/^/2. Then, for each i > 1, the boundary condition Bl is 

= 0, (128) 

r—R 

where Va denotes the covariant derivative with respect to the four metric and where, in terms of the Weyl tensor 
Cabcd, the Newman- Penrose scalar \I'o is given by ^Pq = Cabcdl°"m^l'^'m'^ ■ For L = 1, this reduces to the freezing 
boundary condition proposed in [l^, [13; ll^ i22<l ■ The new boundary conditions (|128|l are local in time and space, 
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and do not depend on the spherical harmonic decomposition. Although they require higher order derivatives of the 
fields at the boundary, high-order derivatives can be eliminated by introducing auxiliary variables at the boundary 
(for example, see Ref. 0). 

Additionally, we have calculated reflection coefficients which quantify the amount of spurious radiation reflected into 
the computational domain both by our new hierarchy of boundary conditions Bl together with constraint-preserving 
boundary conditions (CBPC), and by CPBC currently in use, which freeze the Newman-Penrose scalar to its initial 
value. Including corrections for backscatter, our new boundary conditions, although no longer perfectly absorbing, 
give a reflection coefficient for odd-parity quadrupolar radiation which is less than the one for the freezing condition 
by a factor of M/R for kR > 1.04. (We expect a similar result to hold for even-parity quadrupolar radiation.) 

An application of our results to simulations of the full nonlinear Einstein equations requires that: (i) the spacetime 
near the outer boundary of the computational domain be accurately described by the linearized field equations, (ii) 
the cross sections of the outer boundary surface with the foliation St be approximate metric two-spheres of constant 
area, (iii) the foliation Et near the outer boundary resemble the t — const, foliation of Minkowksi space, where t is the 
standard Minkowski time coordinate, and (iv) the magnitude of the normal component of the shift vector at the outer 
boundary be small compared to one. Criteria (i) and (ii) are fully justified because modern numerical relativity codes 
can push the outer boundary into the weak field regime by using mesh refinement, and can handle spherical outer 
boundaries using multi-block finite differencing [3^l40ll4]| or pseudo-spectral methods [Uli^l. However, criteria (iii) 
and (iv) arc more restrictive because they place requirements on the coordinate and slicing conditions. For example, 
using maximal slicing or a slicing which insures that the mean curvature rapidly decays to zero as one approaches the 
outer boundary might justify criterion (iii), while forcing the normal component (with respect to the outer boundary) 
of the shift vector to be zero at the outer boundary guarantees (iv). On the other hand, these criteria are not justified 
if hyperboloidal slices are used [3^ l6^ \6§L ^9], where the mean curvature asymptotically approaches a constant, 
nonzero value. It should not be difficult to generalize our analysis to more general foliations of Minkowski spacetime 
using the 2-1-1 split discussed in Sect. Ill Bl 

The new boundary conditions l|128(l constructed in this article should be useful for improving the accuracy of 
binary black hole calculations on finite domains. For example, if the outer boundary is spherical with area AttR^ 
and R > lOOM, then the refiection coefficient for CPBC with freezing is less than 0.1% for quadrupolar waves 
with wavelength lOOM or smaller. Since the energy fiux scales as the amplitude of the wave squared, this reflected 
false radiation causes a relative error in the energy flux calculation for quadrupolar radiation of the order 10~® or 
less. If one uses the improved boundary condition B2 proposed in this article instead of the freezing condition, 
then the reflection coefficient is 100 x smaller; ie. less than 0.001%. Correspondingly, the contribution of reflected 
artefactual radiation to the relative error in the energy flux calculation for odd-parity quadrupolar radiation is below 
10""'^'^. Finally, the improved boundary conditions presented here may be useful for minimizing reflections of "junk" 
radiation present in the initial data. 

We would like to conclude by emphasizing two points. The first point is that the hierarchy of new boundary 
conditions Bl proposed in this article (|128|l are not restricted to the Bianchi equations, but can be applied to any 
formulation of Einstein's field equations. It is important that they are used in conjunction with CPBC and suitable 
boundary conditions which control part of the geometry of the outer boundary surface (in particular, they should 
insure that the outer boundary remains spherical and that its area does not change too much in time). While these 
last two types of boundary conditions depend explicitly on the formulation, condition l|128() does not. After all these 
boundary conditions have been specified, one still needs to show that the resulting IBVP is well posed. This issue is one 
which we will address elsewhere. The estimates of the reflection coefficients for spurious gravitational radiation given 
in this article are valid for any representation of the Einstein equations which implements the freezin g \l/n boundar y 
condition together with CPBC. In particular, they are directly applicable to the formulations in lla IitI llSl Il9l l21| . 

The second point is that our improved boundary conditions may not be transparent enough to model accurately all 
physically interesting scenarios on an unbounded domain. For example, it is likely that even with our new boundary 
conditions, one will find an incorrect tail decay when measuring the decay of solutions at a fixed location near the 
outer boundary. The failure of the simple Sommerfeld condition to correctly simulate tail decays for a spheric ally 
symmetric scalar field about a Schwarzschild black hole was demonstrated numerically in 70] . In fact, the work in |7ll | 
proves analytically that the boundary conditions in [70| lead to decay which is faster than any power of 1 /t (whereas 
the expected rate of decay is 1/t^ UM)- future work, we plan to explore ways to improve our new boundary 
conditions to reproduce correctly the tail decay. One possibility is to use the work in ,8, 9', IQj to construct boundary 
conditions which are perfectly absorbing when backscatter is considered. 
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APPENDIX A: STABILITY OF THE ABSORBING BOUNDARY CONDITIONS 

In this appendix, we consider the initial-boundary value problem 



dt-di 



(t){t, r) = 0, t > 0, RQ<r <R, (Al) 



(/.(0,r) =/(r), 9*0(0, r) =g(r), Ro < r < R, (A2) 
(6+)^+i9t0(i,i?o) = 0, =0, i>0, (A3) 

where £, L are natural numbers, < Rq < R are the inner and outer radii of a spherical shell, / and g are smooth 
initial data, and b± — r'^(dt^dr). The reason for introducing the inner boundary at r = i?o is to excise the coordinate 
singularity at r = 0. This is not a restriction for the purpose of this article, since we are interested only in the region 
near the outer boundary, and since we consider the linearized equations for modeling a physical scenario away from 
the strong field region. 

In order to show that the problem (|A1IA2IA3|I is stable in the sense that the solution depends continuously on the 
data, we introduce the following notation: 

$f) = (6±)^V,, J ^0,1,2,3,... 

Notice that $0^'' — ^0 ^ — 4'- By applying the operators &+ and &_ to both sides of Eq. ljAl|l . one finds the formula 

6^(<i>f = ±2jr<i>f ) -{i-j + l){£ + j)rH\^\ 

for j = 1, 2, 3, .... Using this and b±{<^f ^) = <^'-^{, we find that 

2r2at($f ^) = ib+ + 6_)$f ' = ±2jr$f ' - {£ - j + !){£ + j)rH'^^\ + <i>''^\ 

for all j = 1,2,3, .... Therefore, the evolution equations I|A1|I and the boundary conditions (|A3|) yield the evolution 
system 

1 ^-.(+)^^(-)" 



O^'^l^' - (A5) 

- i^a-^-'-^^^^^^^n^ (A6) 

(a, + a,)$i++\ = ^(£±ll$(+\_(£_L)(^ + i + i)$w, (A7) 

(a* - = -^(^-tll<i>W _ (^ „ L){£ + L+ l)4-\ (A8) 



with boundary conditions 



9t$(+\(i,i?o) =0, 9t$^(t,i?) = 0, t>0. (A9) 



The system (|A4IA5IA6IA7IA8|) constitutes a symmetric hyperbolic system with maximally dissipative boundary con- 
ditions IjAOp . It is well-known (see, for example, Ref. 73]) that such systems are well posed and admit energy 
estimates. For example, it follows that a smooth enough solution satisfies the estimate 

E{t) < e^^E{Q) (AlO) 
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with the energy norm 




.2(LH 



L+l 



+ ($(-) (i^^))2 



where 5 is a constant that does not depend on the sohition. In particular, the inequahty (|A10p imphes that the 
solutions depend uniquely and continuously on the initial data. The existence of solutions (including for evolution 
equations with more general potentials than the one in Eq. (lAlfl ) can be proved using methods from semigroup 
theory; see for example chapter 6.3 in [t^ for a well posedness proof for a similar problem. A different well posedness 
proof based on the verification of the Kreiss condition is given in Q . 

Since for £ < L, the exact outgoing solutions cf> ^^g{t,r) constructed in Sect. IIII Bl satisfv the boundary conditions 
(jA3|) . provided the function Ui is compactly supported in {Rq,R), it follows that the boundary conditions (|A3p are 
perfectly absorbing for all £ < L. 
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